Astrophysical Bulletin, vol. 62, No.l, 2007, pp.3-16 February 5, 2008 

Translated from Astrofizichcskij ByuUeten, vol.62, No.l, 2007, pp. 5-19 



Formation of lonization-Cone Structures in Active Galactic 
Nuclei: I. Stationary Model and Linear Stability Analysis 

V. L. Afanasiev\ S. N. Dodonov^ S. S. Khrapov^, V. V. Mustsevoi^, A. V. Moiseev^ 



o 
o 

(N 

cn 



> 

cn 
\o 

(N 
O 

o 

Oh; 

6 



^ Special Astrophysical Observatory, RAS, Nizhnii Arkhyz, Karachai-Cherkessian Republic, 357147 Russia 
^ Volgograd State University, Volgograd, 400062 Russia 

received:September 21, 2006/revised: November 24, 2006 

Abstract. We discuss causes of the formation of the observed kinematics and morphology of cones of ionized 
matter in the neighborhood of the nuclei of Seyfert galaxies. The results of linear stability analysis of an optically 
thin conic jet where radiation cooling and gravity play an important part are reported. The allowance for radiation 
cooling is shown to result in strong damping of all acoustic modes and to have insignificant effect on unstable 
surface Kelvin-Helmholtz modes. In the case of waveguide-resonance internal gravity modes radiative cooling 
suppresses completely the instability of waves propagating away from the ejection source and, vice versa, reduces 
substantially the growth time scale of unstable sourceward propagating modes. The results obtained can be used 
to study ionization cones in Seyfert galaxies with radio jets. In particular, our analysis shows that surface Kelvin- 
Helmholtz modes and volume harmonics are capable of producing regular features observed in optical emission-line 
images of such galaxies. 



1. Introduction 

Observations of many Seyfert galaxies in emission lines 
{Ha, [OIII], etc.) suggest that Narrow Line Regions 
(NLRs) in these galaxies often have conic shapes. 
The [OIII]/Hq line intensity ratio maps constructed by 
some authors (see, e.g., |Pogge, 1989[ |Falcke et al., 1998] ) 
also suggest that gas in NLRs is located inside cone- 
shaped (sometimes bipolar) structures with opening an- 
gles 30° — 110° and linear sizes ranging from several 
tens of parsecs to 20kpc (see Table 3 in the paper of 
[Wilson fc Tsvetanov, 1994| . 

A number of authors thoroughly analyzed the 
structure of NLR in individual galaxies, such as 
Mrk 3 dPogge k De Robertic, 19931 [Capetti et al., 1998| , 



Mrk 573 dTsvetanov fc Walsh, 199"2| [Ferruit et al., 1999 



NGC 3516 dMiyaji et al., 1992[ fVeilleux et al., 1993 



NGC 5252 dWilson fc Tsvetanov, 1994| Moorse et al., 
1998 ), and ESQ 428-G14 dFalcke et al., 19961 ). The 
above authors point out good correlation between 
the directions of ionized cones and those of ra- 
dio jets — their symmetry axes coincide to within 
5 - 10° dWilson fc Tsvetanov, 19941 jFalcke et al., 1996 
jNagar et al., 1999[ ). Figure [T] shows, as an example, two 
[OIII] emission-line images of ionization cones in the 
galaxies NGC 3516 (Sy 1 type) and NGC 5252 (Sy 1.9 
type) and the orientation of the axes of their radio jets. 
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The orientation of the cones with respect to galaxy 
disks differs from one galaxy to another, however, Wilson 
& Tsvetanov (1994) point out that in late-type galaxies 
the cone symmetry axes are virtually perpendicular to the 
disk plane, whereas in early-type galaxies they are inclined 
at small angles to the disk. 

The profile of narrow emission lines in the central re- 
gions of the cone has a complex multicomponent struc- 
ture. This is indicative of several systems of gaseous clouds 
being observed at the corresponding locations along the 
line of sight with mutual velocities as high as several 
hundred fcms"^ ([Capetti et al., 1998} jKaiser et al., 2000} 
jEmsellem et al., 2006| ). In a number of cases, emission- 
line images of the galaxies considered exhibit a Z(S)- 
shaped pattern (Fig. [1]) or other regular structures. 
High spatial resolution HST images of Seyfert galaxies 
dCapetti et al., 19961[Falcke et al., 1998D indicate that the 
Z-shaped pattern begins in the circumnuclear region (at 
galactocentric distances of 10-100 pc) and extends out to 
much larger galactocentric distances (several kpc). The 
NGC 3516 galaxy (Fig. (Ha, and jMiyaji et al., 1992[ ) shows 
Z-shaped filaments in the region r < 10 — 14" (~ 2 kpc) 
with the Northeastern part of the cone extending out to 
more than 6 kpc. The NLR in the NGC 5252 galaxy 
(Fig. [T|d and jMorse et al., 1998^ also has a symmetric Z- 
like shape in the central region r < 10 — 15" (5 — 7) kpc, 
whereas the biconical structure itself, which contains reg- 
ular emission-line "arcs" can be observed out to 18 kpc 
from the center. 
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Fig. 1. The [OIII]-linc images of ionization cones with Z-shaped patterns obtained from observations carried out at 
the 6-m telescope of the Special Astrophysical Observatory of the Russian Academy of Sciences ( [Moiseev et al., 2000[ ). 
The line indicates the orientation of radio jet according to Miyaji et al. (1992) and Wilson & Tsvetanov (1994). (a) 
NGC 3516, (b) NGC 5252. 



The motions of ionized gas inside cones are com- 
plex and poorly studied. Two-dimensional velocity fields 
have been constructed in a few cases ( Veilleux et al., 1993| 



Ferruit et al., 1999 Moiseev et al., 2000), which show 



that different portions of Z-shaped patterns can be both 
"blue-" and "redshifted" with respect to the nucleus. 
Moreover, in some galaxies, like, e.g., in NGC 5252, the 
velocity fields in the region of "arcs" at large distances 
from the nucleus exhibit emission-line filaments that are 
blucshifted exclusively in one cone and redshifted exclu- 
sively in the diametrically opposite cone. 

Ionization cones are believed to be due to the coUima- 
tion of ionizing radiation by the torus of matter accret- 
ing onto a supermassive black hole at the nucleus of the 
galaxy. However, such a scenario of cone formation fails to 
explain the presence of regular structures. Wilson (1993) 
argues that ionized matter moves away from the galac- 
tic center, i.e., that it constitutes a weakly coUimated jet. 
Capetti et al. (1996) carried out a detailed spectrophoto- 
metric study of such a Z-shaped emission-line region in the 
Mrk 573 galaxy and pointed out that radiation of the nu- 
cleus is evidently insufficient to produce the observed NLR 
and that an additional local ionization source is required. 
Ferruit et al. (1999), who used panoramic spectroscopy 
to analyze this object, also concluded that in Mrk 573 the 
necessary additional contribution to ionization is provided 
by shocks produced by the intrusion of the jet from the 
active nucleus into the surrounding clouds of interstellar 
gas. 

There is no consensus of opinion as to what causes 
the formation of regular structures observed in ioniza- 
tion cones. The model of the interaction of the jet with 
gaseous clouds in the circumnuclear region developed by 
Rossi et al. (2000) explains a number of morphological fea- 



tures, but fails to describe the development of symmetric 
Z-shaped features. 

A number of researchers (see, e.g. |Veilleux et al., 1993[ 



Steffen, 1997) believe such features to be hehcal shocks, 



which are due to the presence of a highly coUimated thin 
precessing jet. However, the hypothesis that this is the 
case for all the objects discussed here appears to be too 
daring to say at least. Mulchaey et al. (1992) interpreted 
the Z-shaped structure in NGC 3516 in terms of the bipo- 
lar outflow model where ejected gas is deflected toward 
the galactic disk. Morse et al. (1998) explain the kine- 
matic structure observed in NGC 5252 by the presence of 
three ionized gas disks rotating in differently tilted planes. 
In addition to the above scenarios there is the possibility 
that excitation of a helical shock may be due to shear in- 
stability, which develops at the layer between the matter 
of coUimated radio jet and ionized gas moving at different 
velocities ( [Falcke et al, 1996] ). 

We begin a series of papers where we naturally show 
with no additional assumptions that the observed struc- 
tures and velocity fields can be explained in terms of 
the following scenario, which agrees with the unified 
model of the activity of galactic nuclei ( [Antonucci, 1993[ 
[Wills, 1999P : 



— the highly coUimated high-velocity bipolar jet (ra- 
dio jet) breaks through the torus of optically opaque 
matter accreting onto the supermassive central object 
(black hole) in two diametrically opposite directions 
parallel to the proper angular momentum of the torus 
matter; 

— the jet matter squeezed by the ambient pressure heats 
up intensively and expands rapidly toward the initial 
ejection through the narrow channel thus produced; 
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— the internal gravity waves propagating at an angle to 
the jet axis undergo resonance superreflection at the 
velocity shear surface made up of the jet boundaries; 

— the harmonics of internal gravity waves resonate be- 
tween the jet boundaries and propagate inside the jet 
like in a waveguide; during this process, the energy of 
gravity waves increases with time, amplified via super- 
reflection (resonance waveguide instability); 

~ the development of instability results in the formation 
of a system of nonlinear waves around the jet, which 
heat up the ambient medium; it is important that the 
wave resistance of the ambient medium is significantly 
higher than that of the jet matter and therefore the 
heating mentioned above occurs inside a cone with a 
limited opening angle (which depends on particular 
parameters of the system) around the jet; 

— the allowance for the possible nonlinear superposition 
of different modes and projection effects permits ob- 
taining the qualitative pattern of the observed mor- 
phology of real objects and of the velocity field inside 
the ionization cones. 

As far as we know, no one has yet explored the 
possibility of the development of the above modes, un- 
like the unstable acoustic modes of jets emerging from 
young stellar objects, which were studied by a number 
of authors (see [Ferrari et al., 1982 Payne & Cohn, 1985 
Hardee & Norman, 1988; [Norman fc Stone, 1997| Norman 



& Hardee, 1988). In this paper, we show the principal pos- 
sibility of the buildup of both helical and pinch waveguide- 
resonance internal gravity modes in conic jets, and discuss 
the results of nonlinear numerical modeling of this process. 

In Section [2] we describe the equilibrium model em- 
ployed; in Section [3] we give the linearized equations and 
formulate the problem of determining the eigenfrcqucn- 
cies of unstable jet modes; in Section [3] we discuss the 
dispersion of perturbations during the linear stage of in- 
stability, and in Section [5] we summarize the main conclu- 
sions and make the final comments. We analyze the re- 
sults of nonlinear 2D- and 3D- modeling in our next paper 
( [Afanasiev et al., 2007[ hereafter referred to as Paper II). 

2. Equilibrium model 

Analyses of the dynamics of jet outflows from active galac- 
tic nuclei naturally reveal three characteristic regions in 
radial coordinate r: 

— at r < (1 — 10) pc the gravitational field is determined 
mostly by the Newtonian potential of the central mas- 
sive object; 

- at (1 - 10) < r < (500 - 1000) pc the jet is immersed 
in dispersed mass of the stellar bulge, which can be 
considered to be spheroidal to a first approximation. 
This region coincides with the rigid-rotation region of 



the galactic disk ( jSofue fc Rubin, 2001 ) and therefore 
the gravitational potential can be assumed to be pro- 
portional to squared radius to a fair accuracy; 



— at r > (0.5 — 1) kpc the variation of the gravitational 
potential along the jet outflow differs substantially for 
different galaxies. Moreover, it also strongly depends 
on the orientation of the jet relative to the symmetry 
plane of the galactic disk. 

In this paper we analyze the spectrum of unstable 
modes of a jet located in the gravitational field with poten- 



tial ^ 



Note that wave structures that form in the jet 



in the inner part of the galaxy, where potential can be fit- 
ted fairly well by the Newtonian formula, are incapable of 
distorting the wave pattern over the jet region considered. 
First, because the radial wavelength of perturbations in 
the field of a point mass must decrease with distance from 
the gravitating center, oc r~^/^ ( [Levin et al., 1999| , 
and the spatial scale length is incomparable with the char- 
acteristic wavelength. Second, according to Levin et al. 
(1999), the jet region, where energy can be fed to per- 
turbations via resonance superreflection and hence where 
wave structures with appreciable amplitudes can exist, has 
a limited extent. 

Our aim is to find out whether it is in principle pos- 
sible for unstable modes to develop at the layer between 
the jet and the ambient medium, and therefore we do not 
incorporate the galactic disk into our equilibrium model. 
We show in our next paper (Paper II) that in the case 
of ^(r) oc perturbations in the ambient medium are 
localized in the conic domain near the jet. Hence our for- 
mulation of the problem is formally correct, at least for 
Seyfert galaxies with powerful bulges and ionization cones 
lying outside the plane of the galactic disk. 

We perform our analysis in the spherical coordinate 
system {r,9,ip), where the 9 = axis coincides with the 
symmetry axis of the jet with an opening angle of 9j and 
outflow velocity of V = Vje^. Here is the unit basis vec- 
tor. We model the medium by ideal gas with the following 
equation of state 



(1) 



where pi and pi are the unperturbed (equilibrium) pres- 
sure and density, respectively; Ci is the adiabatic sound 
speed; subscript i is equal to "j" and "a" inside and out- 
side the jet, respectively; we assume that adiabatic index 
7 is the same for the matter of the jet and the ambient 
gas. We assume that gravitational field is spherically sym- 
metric with the center located at the coordinate origin. 
Variations of the gravitational potential can be written in 
the following form: 



(2) 



where D, = const is the angular velocity of gas rotation in 
the circumnuclear region of the disk and = const is a 
normalizing constant. 

We assume that gas outside the jet is at rest. We take 
into account the possible heating of gas of the jet by the 
radiation of the nucleus: qj > 0, where qj = F — pjA 
is the amount of energy absorbed by unit mass per unit 
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time; F = r(r) and A = A(r) are the heating and coohng 
functions, respectively, which depend on temperature T 
exclusively. In the equilibrium state (7^ = outside the 
jet. 

Thus the spatial distribution of the model parameters 
that characterize unperturbed flow has the following form: 



V;p; c;q 



;p„(r);c„(r); 0, e>6 



(3) 



We assume that the jet is contained by the pressure 
of ambient gas and hence that the following equality is 
satisfied at 6* = 



^3- 



PjWcjW = Pa{r)cl{r). 



(4) 



Note that below we make virtually no direct use of 
relation its fulfillment is required for realization of a 
flow with Vg = 0. 

Radial dependences in formula (3) are determined first 
and foremost by the unperturbed balance of forces. Under 
the adopted assumptions, the r component of Euler's 
equations implies 



2 dr 



1 dpt 
Pi dr ■ 



(5) 



It follows from the continuity equation that 

PjVj = AjA^- (6) 

In equation ^ ~ const is the mass-loss rate by the 
system into a unit steradian. It is a free parameter of our 
model. 

Finally, in the case considered the equation of energy 
balance, with formulas ^ and ([6]) taken into account, can 
be written in the following form: 



7-1 



(7) 



Equation of state (P) closes equation set d!])-®. 

We seek the solutions of this equation set in the power- 
law form: /(r) oc r"-'' , where / is any of the parameters 
that characterize the system. It follows from equation ([5]) , 
with equation ((!]) taken into account, that ayV^ -l-VL^r'^ = 
—apC^I^. Hence we find: 



ay = Q^c = 1, oip = ^3, Up = —1, a„ ~ 2. 



(8) 



In this case, for a spherically symmetric potential the 
velocity of matter in the jet relates to the sound speed in 
the ambient medium as: 



7 ■ 



(9) 



The Mach number of the jet is Af — Vj/cj < 1 , i.e. , the 
jet us subsonic. Note that this condition can be satisfied 
simultaneously with Vj/ca ^ 1, making it theoretically 
possible for shocks due to the development of instability 
at the jet boundary to exist in the medium that surrounds 
the jet, because Ca < cj. 



Given that 7 > 1, our model corresponds to an entropy 
distribution St that is stable against convective motions, 
because 



dSj 
dr 



Pi dr 



lap 



37-1 



> 0. 



(10) 



Substitution of power-law radial dependences into for- 
mulas (O and (O and comparison with formula ([9]) yield: 



7(7 - 1) rqj 
37-1 



r.2 ■ 



(11) 



Hence the velocity of matter in the jet is unambigu- 
ously determined by its temperature and heating due to 
external radiation. Heating is rather important for out- 
flows emerging from active galactic nuclei, because the jet 
is illuminated intensively by the radiation of the nucleus. 

We finally determine, in view of formula ([8]), that the 
following condition must be satisfied for the realization of 
the model constructed above: 



5/2 



(12) 



Here Cr and Ca are constants and is the internal en- 
ergy of gas. For further calculations, it is more conve- 
nient to write heating and cooling in terms of internal 
energy, and not in terms of temperature, by taking ad- 
vantage of the fact that (7 — l)e — RT/ji for ideal gas. 
Note that in the temperature interval T < 10^ the sec- 
ond relation in ([12]) agrees excellently with the dependence 
A(T) oc T^-^'^ that is typical of ionized gas in Seyfert galax- 
ies ( [MacDonald fc Bailey, 1981[ [Norman fc Stone, 19"97| . 
Linear temperature dependence of heating function F in 
formula (jl2p is also a good approximation. 

3. Linearized equations and formulation of the 
boundary-value problem 

Let us now analyze the stability of the model that we 
constructed in the previous section against small pertur- 
bations. We proceed from the following set of hydrody- 
namics equations : 



9V 1 
— + (VV)V = — Vp - V*, 
ot p 

|^ + (VV)p + pdivV = 0, 



'dt 



(VV)£ + (7 - l)e divV = Cr e - Ca pe 



5/2 



(13) 



(14) 



(15) 



We obtain the missing equation to close this set by 
choosing the equation of state in the form p = p{p, S) and 
computing its derivative with respect to time: 



dp _ fdp\ dp fdjA dS_ 
dt ^ \dp) gdt^ \9s) p dt 
2 dp Pi dS 2 



(16) 



dt Cv dt 



c7^ + (7-i)P.^g. 
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Our computations take into account the fact that S ~ 
Cy\a{plp^), Cp ~ Cy = R/fi, 7 = Cp/cv, dS/dt = q/T, 
and use the equation of state of ideal gas in the form 
p = RpT/fi. 

Thus the last equation of the above set acquires the 
following form: 



+ (7 - l)pm {Ct - Caps'/^) 



(17) 



To avoid doubts, note that subscript "z" indicates equi- 
librium stationary parameter values. 

We now use standard procedure of linear analysis to 
substitute pressure, internal energy, density, and velocity 
of the medium in the form f(r, 6, Lp, t) — fi{r)+f{r, 6, ip^ t), 
where |/| <C fi. We assume that conditions ©-((Tl) are 
satisfied to derive the following set of equations describing 
the dynamics of small nonadiabatic perturbations written 
for domains that are homogeneous in 9: 



dVr 


or 




'dt ~ 


or 


dve 




ViVe _ 


dt ~ 


or 


r 


dvp 




^ ViV^ 




or 


r 



1 dp p dpi 
dr 



dr 

J_dp 

Pir do ' 

1 



Pi 



dp 



Pir sind dtp^ 



dp 

dt 



1 d 



— {p,V0 sme) + — (p,v^] 



rsmf 



dp dp 

dl^^'d^- 
, - dpi 

+ 

dr 



- dpi 
dr 

PiCa 



dt dr 

3/2- , 3 1/2 _ 

P+2^i P^^ 



(18) 
(19) 
(20) 

(21) 
(22) 



de yj^^ ~ 

dt ^ dr ^ dr 



--^'-^APtS 
1 



rsmf 



3/2 



(7-1)75 (7 + 1) — 

or r 



s +{-/-!) Si 
d{v0sin9) dv^ 



dvr 
dr 



2Vr 

r 



(23) 



de 



dip 



5/2 ~ 

= -Cae/ p. 



We take equation (jlip into account when deriving for- 
mula from formula 

We seek the solutions for perturbed quantities in the 
following form: 



/(r, 0, if, t) = f{e) r'^f exp {ix{r, t) + im^p}. 



(24) 



Based on the requirement of conservation of the flux of 
perturbation energy through a sphere of arbitrary radius 
(r'^pVr = const) , we find the exponent of the radial depen- 
dences of perturbation amplitudes to be /3/ = a/ — jl. 
Here a/ is, as above, the exponent for equilibrium quanti- 
ties, and a* — 2, the exponent for the gravitational poten- 
tial 12]). It is important that for this value formula ([24]) 



gives exact solutions. Let us now introduce the following 
designations: k = dx/dr, u — —dx/dt. In this case, equa- 
tion set (fT8 |) -([23 |) for the equilibrium model considered 
reduces to the following two ordinary differential equa- 
tions: 



— = a;i(wi -)- bi)pir ^, 



di 

de 



1 



rcui^LUi + (5i) 



K + 



sin^ I 



— - ^cot6l. 

Pr 



(25) 
(26) 



Here coi = cu — kVi is the Doppler shifted perturbation 
frequency and 5i — iVi/r^ ^ is the complex amplitude of 
perturbed Lagrangian ^-displacement of the medium such 
that: 



V0 = — = -i[uj 
dt 



kVi)^ = -iCo,t (27) 
We introduced the following designation in formula (P^ : 



A,,; 



I 
r 

2 x2 



k + 



2i 



Pr{uj'i-5: 
Oj — 2(5. 



+ -(fc + 
r 

iAA (37 



l)(fc 



+ ^) 



Cji - (37 - 



■n, 3/2 
-iLAPtSi 

where 



rpi 

37- 



(28) 



A, 



(7 



1 



7r2 C)i + 5i 



-1 



5^) 



3 3/2 



UJi 



3 3/2 ■ 



To solve the boundary-value problem for the determi- 
nation of unstable-mode eigenfrequencies, four boundary 
conditions must be satisfied. We use standard boundary 
conditions at the jet axis based on the physical assump- 
tions of Hardee & Norman (1988) and Norman & Hardee 
(1988). All displacements for axisymmetric perturbations 
should be equal to zero at the jet axis, ^(0) = 0. In the 
case of nonaxisymmetric perturbations (m > 0) the term 
with p(0) in formula (j26p is finite along the jet axis only 
for p(0) = 0. Hence: 



(29) 



m 


= 0, 


m 


> 1, 


m 


-0, 


m 


= 0, 


p{tt) 


= 0, 


m 


> 1, 


e(vr) 


= 0, 


m 


= 0, 



(30) 



To find the boundary conditions at the external {6 = 
9j) jet boundary, we integrate equations (^5)1 and ((^ 
from 9j — e to dj + e (where e ^ 0). Here we adopt 



Pi{0) = Pj - (Pj - Pa)0(9 ~ 9j), where 



9j) is the 



symmetric Heaviside step function. In the same way we 
determine Vi{9), pi{9), ei{9), and Ci{9). Integration yields 
the following conditions: 



p{9,-0)^pi9,+0), 



(31) 
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because tlie riglit-liand sides in equations and (pS)) 
contain no 5-functions (i.e., derivatives of ^-functions) . 
From the physical viewpoint this means the continuity 
of total pressure (j>i + p) at the jet boundary bended by 
perturbations (with the allowance for equation (|4])) and 
the absence of cross-boundary gas flows. 

The dispersion properties of small perturbations in the 
system considered are fully characterized by the following 
dimensionless parameters: 

— constant Mach number M = Vj/cj along the jet; 

— density difference between the ambient gas and jet 
matter p^/p^ ^cpcl; 

— radius-independent fractional radial wavenumber fcr; 

— helical-mode number (the number of spiral arms in az- 
imuth) m; 

— the ratio of the acoustic wave period to the relaxation 
time scale of the gas of the jet r = C\pje'^J^ /kcj] 

— dimensionless phase velocity of perturbations along the 
jet axis z — uj/kcj, which is a solution of equation set 
(l25])-(l26l) with boundary conditions ((291) -([Sll). In this 
case, z is also radius independent. 

The condition of growth of fractional perturbation am- 
plitude {Im z > 0) means that the mode considered is 
unstable. 

Only one of the two parameters — M and R — is inde- 
pendent, whereas the second parameter can be computed 
by relation ^ : = (1 — 1 /R) / 7. We further assume, for 
the sake of simplicity, that gas has the same composition 
throughout the entire system, implying that the external 
cooling coefficient is related to the cooling coefficient 
in the jet as r = y/Rxa- 

Note that in the case of fixed medium with no radiative 
cooling {Vi = 0, Ca = 0, Si = 0, tD, = uj) we set vg = 0, 
m = to immediately derive from the linearized equation 
set for perturbations with wave vector k || Gr the follow- 
ing dispersion equation, which is equivalent to Ai = in 
equation (p6)l . Its solution has the form: 



± kc 



1 



47 - 1 



(33) 



jk^r^ 

Thus in the short- wavelength approximation (fcr ^ 1) 
we obtain the common dispersion law for acoustic waves: 
oj ~ ifccj. The opposite limiting case fcr <C 1 leads to the 
following dispersion relation: 



OJ 



± 



47 - 1 Ci 
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(34) 



Relation (|34|) describes the dispersion law for long- 
wavelength gravity waves. 

4. Discussion of the results of linear analysis 

4.1. Medium Without Relaxation 

To make it easier to identify the causes of the develop- 
ment of unstable modes, we first consider the case of 



a medium without relaxation: r = (adiabatic pertur- 
bations). We use the shooting method to solve numeri- 
cally the boundary-value problem formulated in Chapter 
m Figs.[2]H]show the resulting dispersion curves. The spec- 
trum of unstable modes is discrete and rather complex. 

First, the surface modes resulting from the develop- 
ment of the Kelvin-Helmholtz instability (KHI) in the 
domain of tangential velocity discontinuity between the 
jet and the ambient medium. These modes decay expo- 
nentially with the distance from the jet boundary on both 
sides in the ^-coordinate. 

Second, the modes of the waveguide made up of the 
jet boundary, which are characterized by two "quantum" 
numbers — rij and m. Here Uj is the number of nodes 
of the eigenfunctions of perturbed pressure between the 
boundary of the jet and its symmetry axis determined 
in the direction perpendicular to this axis and m is the 
number of zero points along the jet azimuth (the number 
of arms of the helical spiral in the cross section of the jet). 
Axisymmetric modes with m — are called pinch modes 
and nonaxisymmetric modes are called m-helical modes 
according to the number of zero points in azimuth. The 
harmonics with no zero points along the jet radius are 
the main harmonics and the remaining ones are reflective 
harmonics. 

Third, the spectrum of modes contains weakly unsta- 
ble or decaying ( Imuj < 0) harmonics. For these modes 
a change of any parameter has no effect on the number of 
zero points (i.e., the Ua of this function in the medium sur- 
rounding the jet) of the distribution of perturbed pres- 
sure between the jet boundary and the 9 = tt axis, whereas 
the number Uj of zero points for the jet does change. Such 
modes are a direct consequence of the idealized formula- 
tion of the problem, because the ambient medium can also 
be formally viewed as a waveguide. However, in the real 
situation the development of a standing wave in 6 coordi- 
nate between the jet boundary and the 9 = tt axis shall 
be impeded by local heterogeneities present. 

In the situation considered the dispersion properties of 
perturbations differ substantially from those in the case 
of the gravitational potential of a compact object. This 
is due, first and foremost, to weak compressibility of the 
medium in the jet on the one hand, and to appreciable 
density stratification on the other hand. 

The local dispersion law allows the existence of two 
types of oscillatory modes in each of these media. These 
are gravity acoustical waves (GAW) and internal gravity 
waves (IGW). GAWs are common longitudinal acoustic 
waves modified by gradient effects. IGWs are due to the 
shear elasticity of the medium resulting from the disbal- 
ance of buoyancy and gravity forces. In the limiting case of 
incompressible medium {M <^ 1) they become transverse 
waves. 

Our model has a preferred direction — that of the ve- 
locity shear vector, which is parallel to the gravity force, 
— and therefore for each of these types of modes one can 
single out waves propagating in the direction of the ve- 
locity shear and in the opposite direction. Therefore, as is 
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Fig. 2. Dimensionless phase velocity Re{u} /kcj) (a, c) and amplitude increments Im (uj/kcj) (b, d) as functions of 
the Mach number M = Vj/cj for axisymmetric (a, b) and helical m ~ 1 unstable modes (c, d). The number nj of zero 
points of eigenfunctions inside the jet is indicated near each curve. The jet half-opening angle is equal to dj — 20°, 
kr = 5. The dashed and solid lines show the and u~ mode families, respectively. 



evident from Fig. [2KjC, all modes break into two families: 
in the reference frame connected with matter modes of 
one family (w^) propagate away from the ejection source 
and modes of the other family (u^) propagate toward the 
source. In the case of a medium with no radiative cooling 
all these modes turn out to be unstable over a wide range 
of parameters. 



The short-wavelength approximation of the corresponding 
dispersion equation has the following form (see Appendix): 



-gkj 



din po 



dz 



{kl + kl)cl 



, In Po 

dz 



= 0, 



UJ 



(35) 



One can determine that most of the dispersion curve 
shown in Fig.[2]corresponds to IGW by directly comparing 
our results with the solution of the well-known problem of 
wave dispersion in a compressible medium with vertical 
density stratification due to uniform gravity field —gez- 



kl. 



where a) = w — kV, k^ = k^ 

We now use the short-wavelength approximation to 
make the formal transition from equation (j35p to the case 
considered by substituting the r coordinate for the z coor- 
dinate, k for kz, k\ = k"^ -\- m? /r^ for fc^, pi for po, Q for 
Cs, and Pi for po- Given the power-law behavior of the de- 
pendences of thermodynamic parameters with exponents 
dH), we normahze equation ([35]) by k'^cf to find the local 
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Fig. 3. Dispersion curves from Fig. [2^ (solid curves), 
gravity-acoustic (GAW) and internal gravity (IGW) eigen- 
modes of the jet and atmosphere determined from eq. (|36p . 
([57|) (dashed curves). 

dimensionless dispersion laws for the ambient medium and 
jet, respectively: 



{z - Mf 



I + 5^ + 



{z - Mf 



■(37-l)-0, 



(37 



- 1) = 



{Rzf 



(36) 



(37) 



Parameter d = k±/k characterizes the inclination of the 
perturbation propagation vector with respect to the radial 
direction. 

Figure [3] shows all solutions of equations ([5^ - ([57)) . the 
straight line z — M corresponding to the velocity of the 
jet matter, and the dispersion curves from Fig. [2^. It is 
evident that the dispersion curves of our jet modes tend 
to the curves corresponding to IGWs. 

Like in the case of jets emerging from young 
stellar objects ( [Ferrari et al., 1982{ [Payne fc Co hn, 1985*; 
[Hardee fc Norman, 1988[ [Norman fc Hardee, 1988) , in- 
stability in our case is due to superreflection and super- 



refraction (Miles, 1957 Ribner, 1957) 



Flow in the jet is subsonic in our model. However, the 
excess of velocity shear at the jet boundary over the wave 
velocity along this boundary, which is required for super- 
reflection, is achieved owing to the smallness of the char- 
acteristic propagation velocity of internal gravity waves. 
For the parameter values corresponding to the curves 
shown in Figsj^HH this corresponds to an ~ 0.2 increase 



in the Mach number. Thus the allowance for gravity re- 
sults in the appearance of additional unstable jet modes 
— waveguide-resonance internal gravity waves whose am- 
plification mechanism is due to the superreflection of this 
type of waves from the jet boundary. 

Because of their low Mach numbers, gravity-acoustic 
waves do not satisfy the superreflection condition. In this 
case the main physical cause of instability is the Bernoulli 
effect — the well-studied Kelvin-Helmholtz instability de- 
velops at the velocity-shear layer between the jet matter 
and the surrounding gas. At the same time, these waves 
resonate between the jet boundary and the 9 = n axis and 
the process is accompanied by the formation of a weakly 
unstable wave in the case of a medium without relaxation. 

The amplitude increment of each reflective harmonic 
{rij > 1) grows rapidly with decreasing radial wavenum- 
ber (see Fig.|4|). However, this is not a physical effect, but 
only a result of normalization by kcj. It becomes clear 
below (see 14. 2|) that renormalization to the wavenumber 
independent Vaisala frequency results in the unlimited de- 
crease of the growth time scale of unstable IGW modes 
with decreasing radial wavelength. 

An important result is that the time scales of the de- 
velopment of the instability discussed here depend only 
slightly on the jet opening angle (see Fig. [5|). 

The extremely weak dependence of the increments of 
unstable modes on jet parameters leads us to suggest that 
different modes may develop simultaneously and coexist. 
On the other hand, we can conclude that the wavelengths 
at which the instability builds up, should be to a greater 
degree determined by the initial perturbations. 

At the same time, the amplitude increment of reflective 
harmonics decreases and that of the main harmonics (rij = 
0) — both the axisymmetric and the first helical modes of 
internal gravity waves — on the contrary, increases with 
increasing kr. Moreover, at kr > 10 these modes are most 
likely to develop instability. 

4.2. Effect of Radiative Cooling on the Dispersion of 
Unstable Modes 

The allowance for radiative losses changes radically the 
spectrum of unstable modes in the system considered: 

GAW modes are completely suppressed by radiative 
cooling and their decay decrements exceed the frequency 
by a factor of several tens to several hundred. This actually 
means that such modes cannot even develop. We therefore 
do not show the corresponding dispersion curves in our 
figures in order not to encumber them. 

Surface KHI modes are extremely sensitive to ra- 
diative cooling. Their growth time scale increases only 
slightly compared to the case of adiabatic perturbations 
(see Figs.[i-[7|). 

Waveguide-resonance modes of IGW of the 
u+ family become damped already at small values of 
radiation-cooling parameter r. Their damping time scale 
decreases rapidly with increasing r (see Fig. [7b)- It can 
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Fig. 4. Dependences of dimensionless phase velocity Re{uj/kcj) (a) and fractional growth rate of the amplitude of 
unstable modes Im (uj/kcj) (b) on dimensionless wavenumber kr. The solid and dashed curves show the modes with 
m = and m = 1, respectively. The number of zero points of eigenfunctions inside the jet is indicated near each curve. 
The half-opening angle of the jet is 6^ = 20°, M = 0,7, 




Fig. 5. Dimensionless phase velocities Re{uj/kcj) (a) and amplitude increments Im{Lu/kcj) (b) as functions of the 
jet half-opening angle 9j in degrees. Designations are the same as in Fig. 21 kr — 5, M = 0.7. 



therefore be argued that they are also fully suppressed 
by radiative cooling. We show the dispersion curves cor- 
responding to these modes only in Fig. [T] Note that 
strictly speaking, in the presence of radiative cooling these 
modes should be called entropy-vortex-type and not inter- 
nal gravity modes. We nevertheless retain the old name to 
facilitate the comparison to the case of the medium with- 
out relaxation. 



Waveguide-resonance IGW modes of the u 

family increase in strength substantially because of ra- 
diation losses. Their increment increases with increasing 
T (Fig. [7]) without reaching saturatiorQ. 

To explain the latter effect, we must use the formula 
relating wave energy density E in a moving medium to 



^ We verified this statement by computations up to r = 100 



10 



Afanasiev et al.: Formation of Ionization- Cone Structures.. I. 



0,7 



0.6 



0.5 



0,2 



0,1 - 



0.0 







(a) 


4^2 

















0.0 



0,3 



0,4 

M 



0,6 



0.8 




Fig. 6. Dimensionless phase velocities Re{uj/kc-j) (a) and amplitude increments Im(uj/kcj) (b) as functions of the 
Mach number for modes with radiative cooling (r = 5). Each curve is labelled by the number of the harmonic (the 
number of zero points of pressure between the boundary and the jet axis). Only u^-family harmonics of mode to = 
are shown. 6*, = 20°, kr = 5. 
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Fig. 7. Dimensionless phase velocities Re{uj/kcj) (a) and amplitude increments Im (uj/kcj) (b) as functions of 
radiative-cooling parameter r for different unstable modes. Each curve is labelled by the number of the harmonic 
(the number of zero points of pressure between the boundary and the jet axis). The dashed lines show the modes of 
the family. Only harmonics of mode to = are shown. 9j = 20°, M = 0.5, kr = 5. 



wave energy density Eq in the reference frame comoving Note that although formula (|38p was initially derived for 
with the medium ( [Landau fc Lifshitz, 1987| : acoustic waves, it is of universal nature, because it al- 

lows simple quantum interpretation. Namely, the number 
of photons of the wave field J\f — inE/hu = 2TTEo/h{u! — 
kV) does not depend on the choice of the reference frame 



E — En 



'uj - kV 



(38) dLandau fc Lifshitz, 1987p . 
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Fig. 8. Dependences of dimensionless phase velocities Re{uj/kcj) (a) and amplitude increments Im{Lu/kcj) (b) on 
the jet half-opening angle Oj for different unstable modes. Each curve is labelled by the number of the harmonic (the 
number of pressure zero points between the boundary and axis of the jet). Only harmonics of mode m = are shown. 
kr ~ 5, T = 5. 



Thus modes of the it"*" family, for which Reu> > kVj 
{Rez > M), have positive energy density in the jet, 
whereas modes of the u~ family, for which < Re uj < kVj 
{0 < Re z < M), have negative energy density. Therefore 
the decrease of the wave energy due to radiative cooling 
decreases the energy of u^-family modes, and, correspond- 
ingly, results in their decay. And vice versa, the decrease 
of the wave energy due to radiative cooling increases the 
absolute value of the energy density of u~-family modes 
and results in their amplification. The latter situation is 
a typical example of radiative and dissipative instability. 

Like in the case of a medium without dissipation the 
dispersion law depends only very slightly on the change 
of the jet opening angle over a wide range of its values 
(Fig.®. 

Fig. [9] shows how the frequency of unstable perturba- 
tions depends on dimensionless wavenumber kr. It is con- 
venient to normalize this frequency not to the frequency 
of acoustic waves, which itself depends on fc, but to the 
characteristic frequency of IGWs — the Brunt- Vaisala fre- 
quency: 

^-lfLfl^^_^!^^U^i. (39) 
Pj ar \pj ar pjC^ ar J 7^ 

As is evident from Fig.[5jD, our analysis predicts unlim- 
ited decrease of the characteristic time scale of the growth 
of unstable IGW modes with decreasing radial wavelength. 
One must nevertheless bear in mind that the presence of 
a transitional layer of finite thickness where velocity 
varies smoothly from Vj inside the jet to zero outside the 
jet, stabilizes perturbations with wavelength X < I. 



The characteristic feature of the IGW modes con- 
sidered is that the maximum of perturbed pressure is 
achieved at the jet boundary, whereas both the perturbed 
displacement ^ of this boundary in the direction transver- 
sal to it and density perturbation are equal to zero in the 
adiabatic case and increase insignificantly with increas- 
ing parameter r. Moreover, outside the jet the perturba- 
tion amplitudes decrease very rapidly with the distance 
from the jet due to the difference of the impedances of 
the media: pjCj < PaCa- We show in Paper II that dur- 
ing the linear stage of the development of instability the 
initial flow (i.e., the jet proper) is not destroyed because 
of the smallness of the 6'-displacement, and perturbations 
remain localized inside a cone near the jet. 

In conclusion, we make two small comments. First, al- 
though Figs. [6H9] show the dispersion curves only for ax- 
isymmetric pinch modes, the results for helical (m > 1) 
modes are qualitatively the same. Second, the compu- 
tations for bipolar outflows with the subsequent change 
of the boundary condition ([50]) yield the results that are 
identical to those that we discuss here for unipolar out- 
flow. This conclusion follows directly from the localization 
of unstable modes near the jet boundary that we discuss 
in this paper. 

5. Conclusions 

Our linear analysis leads us to conclude that: 

— Conical mass outflows in a field of quadratic gravi- 
tational potential and similar to those observed in a 
number of Seyfert galaxies are unstable against the 
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Fig. 9. Dependences of dimensionless frequency Reiuj/N) (a) and increment Im{Lo/N) (b) on dimensionless radial 
wavenumber kr for different unstable modes. Each curve is labelled by the number of the harmonic (the number of 
zero points of pressure between the boundary and axis of the jet). Only the harmonics of mode m = are shown. 
M = 0.5, Oj = 20°, r = 5. 



resonance-waveguide development of a wide spectrum 
of pinch and helical internal gravity waves. 

— The characteristic amplitude growth time scale of 
these modes depends extremely slightly on the jet 
opening angle over a wide range of theses angles. 

— Radiative cooling suppresses completely all gravita- 
tional acoustic waves, has only a slight effect on un- 
stable surface Kelvin-Helmholtz modes, and results 
in the decay of waveguide-resonance internal grav- 
ity modes propagating in antisource direction rel- 
ative to the jet matter. And vice versa, radiative 
cooling increases substantially the instability of such 
sourceward-propagating modes. 

— amplification mentioned above has the form of radia- 
tive and dissipative instability of modes with negative 
energy density. 

— the formation of the observed regular patterns in ra- 
diation cones in the vicinity of the nuclei of Seyfert 
galaxies can be due only to unstable surface modes 
and slow (sourceward moving in the jet) waveguide- 
resonance IGWs modes. The velocities of these modes 
along the jet boundary exceed the characteristic sound 
speed in the ambient atmosphere and this fact allows 
us to believe in their possible evolution into shocks. 

— In the case of small jet opening angles the main har- 
monic of the pinch mode of IGWs is most likely to 
develop in the short-wavelength domain (kr > 20). In 
the longer-wavelength domain (5 < fcr < 15) the main 
harmonic of the first helical mode is most likely to de- 
velop. 

— Because of the different localization of the first heli- 
cal and pinch modes the development of one of these 



modes should no fatal consequences whatsoever for the 
other mode. 

We describe detailed numerical simulations of the de- 
velopment of these waves in our forthcoming Paper II. 
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Appendix A: Dispersion relation for IGWs 

Here we briefly describe the derivation of the disper- 
sion relation for internal gravity waves. We consider low- 
amplitude waves in a compressible medium with vertical 
density stratification due to uniform gravitational field 
—ge^, where g = const. The initial linearized set of hy- 
drodynamics equations has the following form: 



dp , . dpo 
at oz 



Po 



div I 



dvz 
dt 



V±p, 

Po 

1 dp p 

— « 6'—. 

Po oz Po 



dp , , dpo 



dt ' dz '"'[dt^'"' dz 



dp , , dpv. 



(A.l) 
(A.2) 
(A.3) 
(A.4) 



where kz is the vertical component of wave vector k (k^ = 
-f k^^). We then seek a solution in the form of planar 
waves: 



/(x, y, z,t) — f exp{ikxX + ikyy + ikzZ — iut}, 



(A.7) 



where / is the perturbed function with amplitude / = 
const. The set of differential equations (jA.l[) - (IA.4|) trans- 
forms into the following algebraic equations: 



- iujp + Vz'^ + ipo(k±v_L + kzVz) = 0, 
dz 



Po 

■ , P P 



Po Po 



tup - pogvz 



. . , . dpo 
lup + Vz —r- 

az 



(A.8) 
(A.9) 
(A.IO) 
(A.ll) 



This equation set is homogeneous and has a nontriv- 
ial solution only if its determinant is identically equal to 
zero. The latter condition yields the following dispersion 
equation: 



~k'cl- 



din pq 
dz 



-g- 



—gk 



jdlnpo 

' dz 



0.(A.12) 



In the case of a medium moving with velocity V the 
Doppler transform for frequency lu — uj — kV yields equa- 
tion (|35p in the main part of the paper. 

Acknowledgements. We are grateful to I. G. Kovalenko for his 
critical comments and to V. V. Levi for numerous useful dis- 
cussions. The images of the NGC 3516 and NGC 5252 galax- 
ies were obtained with the 6-meter telescope of the Special 
Astrophysical Observatory of the Russian Academy of Sciences 
funded by the Ministry of Science of the Russian Federation 
(registration number 01-43). This work was supported in part 
by the "Nonstationary objects in the Universe" program of the 
Ministry of Industry and Science of the Russian Federation. 
A. V. Moiseev and V. L. Afanasiev also acknowledge the of 
from the Russian Foundation for Basic Research (project no. 
06-02-16825). 



Here subscript z denotes vertical components and sub- 
script _L, the vector components orthogonal to e^. The 
equation of hydrostatic equilibrium can be written in the 
following form: 
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1 dpo{z) 
poiz) dz 



(A.5) 



Other designations are the same as in the main part of 
the paper. 

We consider short-wavelength perturbations along the 
z axis: 



1 



^In Po 

dz 



«1, 



(A.6) 



